The aim of this practical work is to get acquainted with mathematical morphology transformations. Applying different operations with several structuring elements on simple images will allow you understanding the actions, effects and properties of the operations, the role of the structuring element, as well as the need for appropriate pre-processing depending on the application at hand.
In this session, we will mainly focus on functions from the skimage.morphology package.
First, you have to download a set of images that you will use throughout the TP: https://partage.imt.fr/index.php/s/9YPTAzKsiXHQWH3
## imports
import matplotlib.pyplot as plt
import skimage.io as skio
import skimage
import numpy as np
import skimage.morphology as morpho
import skimage.segmentation as seg
student = True
The basic idea in mathematical orphology is to probe an image with a simple, pre-defined shape, drawing conclusions on how this shape fits or misses the shapes in the image. This simple "probe" is called the structuring element. For example, it can take the shape of a square, a disk or a diamond, with a predefined radius.
Use the function morpho.disk() of the morpho library to visualize a disk structuring element of radius 5.
if student :
strell = morpho.disk(5)
print(strell)
plt.imshow(strell)
plt.show()
[[0 0 0 0 0 1 0 0 0 0 0] [0 0 1 1 1 1 1 1 1 0 0] [0 1 1 1 1 1 1 1 1 1 0] [0 1 1 1 1 1 1 1 1 1 0] [0 1 1 1 1 1 1 1 1 1 0] [1 1 1 1 1 1 1 1 1 1 1] [0 1 1 1 1 1 1 1 1 1 0] [0 1 1 1 1 1 1 1 1 1 0] [0 1 1 1 1 1 1 1 1 1 0] [0 0 1 1 1 1 1 1 1 0 0] [0 0 0 0 0 1 0 0 0 0 0]]
Find a way to create a vertical line as a structuring element with parameter length l :
def line_strell(width,length):
if student :
strell = np.ones((length, width), dtype=np.uint8)
# Fill the center column with ones to create a vertical line
return strell
strell = line_strell(2,5)
print(strell)
plt.imshow(strell)
plt.show()
[[1 1] [1 1] [1 1] [1 1] [1 1]]
In this section, we will work on a cyto-image. Make sure to properly load the "cell.tif" image. Apply a thresholding of level t = 150 to that image to turn it into a binary image. Visualize the result in greyscale, by using the imshow function.
#load the cell image in grayscale and display it
import cv2
if student:
img = skio.imread("cell.tif")
plt.imshow(img, cmap = 'gray')
plt.show()
_, binary_image = cv2.threshold(img, 150, 255, cv2.THRESH_BINARY)
plt.imshow(binary_image, cmap = 'gray')
plt.show()
We will first start by applying simple binary morphological operations on those images. You will try the following operations:
You will try to visualize the results by varying the structuring element and its size.
strell = morpho.disk(2)
fig,ax = plt.subplots(2, 2, figsize=(10,10))
imgs = [morpho.erosion(binary_image, strell), morpho.dilation(binary_image, strell), morpho.opening(binary_image, strell),morpho.closing(binary_image, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
strell = morpho.disk(5)
fig,ax = plt.subplots(2, 2, figsize=(10,10))
imgs = [morpho.erosion(binary_image, strell), morpho.dilation(binary_image, strell), morpho.opening(binary_image, strell),morpho.closing(binary_image, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
strell = morpho.disk(7)
fig,ax = plt.subplots(2, 2, figsize=(10,10))
imgs = [morpho.erosion(binary_image, strell), morpho.dilation(binary_image, strell), morpho.opening(binary_image, strell),morpho.closing(binary_image, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
# try your own structuring elements and make the size vary
if student :
strell = line_strell(2,5)
fig,ax = plt.subplots(2, 2, figsize=(10,10))
imgs = [morpho.erosion(binary_image, strell), morpho.dilation(binary_image, strell), morpho.opening(binary_image, strell),morpho.closing(binary_image, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
strell = line_strell(4,10)
fig,ax = plt.subplots(2, 2, figsize=(10,10))
imgs = [morpho.erosion(binary_image, strell), morpho.dilation(binary_image, strell), morpho.opening(binary_image, strell),morpho.closing(binary_image, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
...
In this section, we will use the bateau.tif image. Make sure to load it and visualize it in greyscale.
#load the bateau image in grayscale and display it
if student:
img = skio.imread("bat.tif")
plt.imshow(img, cmap = "gray")
plt.show()
Apply the same 4 operations as before (dilation,erosion,opening,closing) and visualize the results. Try out different structuring elements and different sizes.
if student:
### start code
strell = morpho.disk(2)
fig,ax = plt.subplots(2, 2, figsize=(20,20))
imgs = [morpho.erosion(img, strell), morpho.dilation(img, strell), morpho.opening(img, strell),morpho.closing(img, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
strell = morpho.disk(8)
fig,ax = plt.subplots(2, 2, figsize=(20,20))
imgs = [morpho.erosion(img, strell), morpho.dilation(img, strell), morpho.opening(img, strell),morpho.closing(img, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
strell = strell = line_strell(2,5)
fig,ax = plt.subplots(2, 2, figsize=(20,20))
imgs = [morpho.erosion(img, strell), morpho.dilation(img, strell), morpho.opening(img, strell),morpho.closing(img, strell)]# Add in this list the 4 images, obtained with the for operations
names = ["erosion","dilation","opening","closing"]
for i in range(4):
ax[i//2,i%2].set_title(names[i])
ax[i//2,i%2].imshow(imgs[i],cmap = "gray")
plt.show()
### end code
result of a dilation by a disk structuring element of size 2 followed by a dilation by a disk structuring element of size 3? Compare it with a dilation by a disk of size 5. Explain.
, where $\bigoplus$ is the Minkowski sum So in this case the Minkowski sum of a circle with radius 2 with a circle of radius 3, returns exactly a circle of radius 5. So, do this process in sequence should return the same values for the two operations. Problably, we cannot see exactly the same result because the circle is discretized, and not continous.
if student:
### start code
dilation = morpho.dilation(img, morpho.disk(5))
plt.imshow(dilation,cmap = "gray")
plt.title("Single dilation")
plt.show()
sequential_dilation = morpho.dilation(img, morpho.disk(2))
dilation_sequential_dilation = morpho.dilation(sequential_dilation, morpho.disk(3))
plt.imshow(dilation_sequential_dilation,cmap = "gray")
plt.title("Sequential dilation")
plt.show()
plt.imshow(np.abs(dilation-dilation_sequential_dilation),cmap = "gray")
plt.title('Difference between the two images')
plt.colorbar()
plt.show()
print(np.sum(np.abs(dilation-sequential_dilation)))
4966894
an opening by a disk structuring element of size 3? Compare with a dilation of the image by (a disk of size 2 dilated by a disk of size 3). Please note that you might have to use the np.pad function to perform this. Conclude.
Thus, $$D(X, D(B, B')) = D(X, B\bigoplus B') = X \bigoplus B\bigoplus B'$$
And $$D(D(X, B), B') = D((X \bigoplus B), B') = X \bigoplus B\bigoplus B'$$
Same results.
if student:
# First create the sequential structuring element, that is a disk of size 2 dilated by a disk of size 3
sqr_strell = morpho.disk(2)
sqr_strell = np.pad(sqr_strell, 3)
sqr_strell = morpho.dilation(sqr_strell, morpho.disk(3))
#print(sqr_strell)
# Use it to perform the image dilation
dilation = morpho.dilation(img, sqr_strell)
plt.imshow(dilation,cmap = "gray")
plt.title("Single dilation")
plt.show()
# Compare it with the sequential dilation
sequential_dilation = morpho.dilation(morpho.dilation(img, morpho.disk(2)),morpho.disk(3))
plt.imshow(sequential_dilation,cmap = "gray")
plt.title("Sequential dilation")
plt.show()
plt.imshow(dilation-sequential_dilation,cmap = "gray")
plt.title('Difference between the two images')
plt.colorbar()
plt.show()
print(np.sum(np.abs(dilation-sequential_dilation)))
0
if student :
### start code
closed_once = morpho.closing(img, morpho.disk(5))
plt.imshow(closed_once, cmap='gray')
plt.show()
closed_twice = morpho.closing(img, morpho.disk(5))
plt.imshow(closed_twice, cmap='gray')
plt.show()
plt.imshow(closed_twice - closed_once, cmap='gray')
plt.show()
print(np.sum(np.abs(closed_twice - closed_once)))
### end code
0
Comment the result depending on the choice of the structuring element.
So, when we realize a subtraction of the original image and the oppening, we highlight the brighters and thin elements in comparission to the structure element.
if student:
### start code
top_hat = img - morpho.opening(img, morpho.disk(1))
plt.figure(figsize = (7,7))
plt.imshow(top_hat,cmap = "gray")
plt.show()
if student:
### start code
top_hat = img - morpho.opening(img, morpho.disk(2))
plt.figure(figsize = (7,7))
plt.imshow(top_hat,cmap = "gray")
plt.show()
if student:
### start code
top_hat = img - morpho.opening(img, morpho.disk(10))
plt.figure(figsize = (7,7))
plt.imshow(top_hat,cmap = "gray")
plt.show()
Comment the result depending on the choice of the structuring element.
if student:
### start code
strell = morpho.disk(2)
dual_top_hat = morpho.closing(img, strell) - img
plt.figure(figsize = (7,7))
plt.imshow(dual_top_hat,cmap = "gray")
plt.show()
operations?
# Example
img = skio.imread("montreuil.tif")
# Structuring elements
se1=morpho.disk(1)
se2=morpho.disk(2)
se3=morpho.disk(3)
se4=morpho.disk(4)
se5=morpho.disk(5)
# Images
fas1= morpho.closing(morpho.opening(img, se1), se1)
fas2= morpho.closing(morpho.opening(img, se2), se2)
fas3= morpho.closing(morpho.opening(img, se3), se3)
fas4= morpho.closing(morpho.opening(img, se4), se4)
fas5= morpho.closing(morpho.opening(img, se5), se5)
fig,ax = plt.subplots(2,3,figsize = (20,15))
ax[0,0].imshow(fas1,cmap = "gray")
ax[0,1].imshow(fas2,cmap = "gray")
ax[0,2].imshow(fas3,cmap = "gray")
ax[1,0].imshow(fas4,cmap = "gray")
ax[1,1].imshow(fas5,cmap = "gray")
ax[1,2].imshow(img,cmap = "gray")
plt.show()
#Try you own
if student :
se1=morpho.disk(1)
se2=morpho.rectangle(2,5)
se3=morpho.rectangle(3, 6)
se4=morpho.rectangle(4,10)
se5=morpho.rectangle(6,12)
# Images
fas1= morpho.closing(morpho.opening(img, se1), se1)
fas2= morpho.closing(morpho.opening(img, se2), se2)
fas3= morpho.closing(morpho.opening(img, se3), se3)
fas4= morpho.closing(morpho.opening(img, se4), se4)
fas5= morpho.closing(morpho.opening(img, se5), se5)
fig,ax = plt.subplots(2,3,figsize = (20,15))
ax[0,0].imshow(fas1,cmap = "gray")
ax[0,1].imshow(fas2,cmap = "gray")
ax[0,2].imshow(fas3,cmap = "gray")
ax[1,0].imshow(fas4,cmap = "gray")
ax[1,1].imshow(fas5,cmap = "gray")
ax[1,2].imshow(img,cmap = "gray")
plt.show()
dilation and erosion with a structuring element of radius1). What do you observe ?
The result of the morpholical gradient highlights the contours of the cells. This is similar to the norm of the regular gradient.
if student :
img = skio.imread("cell.tif")
erosion = morpho.erosion(img, morpho.disk(1))
dilation = morpho.dilation(img, morpho.disk(1))
morpho_grad = dilation - erosion
plt.figure(figsize=(10,10))
plt.imshow(morpho_grad,cmap = "gray")
plt.show()
As we can see in the previous question the morphological gradient contours even the smaller contours, so it makes that when we threshold it, it captures undesired patterns.
if student :
img = skio.imread("cell.tif")
_,binary_image = cv2.threshold(morpho_grad, 25, 255, cv2.THRESH_BINARY)
plt.figure(figsize=(10,10))
plt.imshow(binary_image,cmap = "gray")
plt.show()
What do you observe ?
The idea of the watershed algorithm is to find out basins of an image, which means, find out pixels that shares the same local minima(based on intensity) in a region. So we can visualize the colors of the regions and their contours.
import matplotlib
from skimage.feature import peak_local_max
cmaprand = matplotlib.colors.ListedColormap ( np.random.rand ( 256,3))
if student :
wat = seg.watershed(morpho_grad, markers =50, watershed_line=True)
plt.imshow(wat, cmap =cmaprand)
plt.show()
morpho_grad_w_edges = np.copy(morpho_grad)
morpho_grad_w_edges[wat == 0] = 0
plt.imshow(morpho_grad-morpho_grad_w_edges, cmap='gray')
plt.show()
It seems that for high values of markers that we start to see some structures that were not merged in the original image and making perceptible the details of some borders.
The result seems worse than the previous operations, but it seems to help in delete fake conexions beetween the cells generate by noise.
wat_closing = seg.watershed(morpho.closing(morpho_grad, morpho.disk(2)), markers = 100, watershed_line=True)
wat_opening = seg.watershed(morpho.opening(morpho_grad, morpho.disk(2)), markers = 100, watershed_line=True)
morpho_grad_opening = np.copy(morpho_grad)
morpho_grad_opening[wat_opening == 0] = 0
morpho_grad_closing = np.copy(morpho_grad)
morpho_grad_closing[wat_closing == 0] = 0
fig, axes = plt.subplots(1,3, figsize=(15,5))
axes[0].imshow(wat, cmap=cmaprand)
axes[0].set_title('watershed')
axes[1].imshow(wat_closing, cmap=cmaprand)
axes[1].set_title('pre-closed watershed')
axes[2].imshow(wat_opening, cmap=cmaprand)
axes[2].set_title('pre-opened watershed')
plt.show()
fig, axes = plt.subplots(1,3, figsize=(15,5))
axes[0].imshow(morpho_grad_w_edges, cmap='gray')
axes[0].set_title('original contours')
axes[1].imshow(morpho_grad_closing, cmap='gray')
axes[1].set_title('pre-closed watershed differnce')
axes[2].imshow(morpho_grad_opening , cmap='gray')
axes[2].set_title('pre-opened watershed differnce')
plt.show()
fig, axes = plt.subplots(1,3, figsize=(15,5))
axes[0].imshow(morpho_grad, cmap='gray')
axes[0].set_title('original contours')
axes[1].imshow(morpho_grad_closing - morpho_grad_w_edges, cmap='gray')
axes[1].set_title('pre-closed watershed differnce')
axes[2].imshow(morpho_grad_opening - morpho_grad_w_edges, cmap='gray')
axes[2].set_title('pre-opened watershed differnce')
plt.show()
Your objective is to try to do better than what you did with the closure and opening. To do that, we propose to use specific markers for the watershed transform.
Q : In the next cell, compute the reconstruction from an opening of the gradient. You can use as structuring element a disk of size 5. What do you observe ? Could it help you to efficiently segment the cells ?
A The image seems to have more noise, but the contours has became more defined. The definition of the countours can help in the segmentation.
strell = morpho.disk(5)
imOpen = morpho.opening(morpho_grad, strell)
reco = morpho.reconstruction(imOpen, morpho_grad)
plt.figure(figsize =(10,10))
plt.imshow(reco,cmap = "gray")
plt.show()
Q : In the next cell, compute the reconstruction from an closing of the gradient. You can use as structuring alement a disk of size 5. What do you observe ? Could it help you to efficiently segment the cells ?
A : The borders are more smooth and in overall the image seems more denoised. Again, it might be usefull in segmentation.
imClos = morpho.closing(morpho_grad, strell)
reco = morpho.reconstruction(morpho_grad, imClos)
plt.figure(figsize =(10,10))
plt.imshow(reco, cmap = "gray")
plt.show()
Q : This last part is less guided. Use the reconstruction operation and the above questions to find good markers for the watershed algorithm. You could benefit from using the peak_local_max function of skimage.feature
from skimage.feature import peak_local_max
plt.imshow(img, cmap='gray')
plt.title('Original')
plt.show()
strell = morpho.disk(1)
img_grad = morpho.dilation(img, strell) - morpho.erosion(img, strell)
plt.imshow(img_grad, cmap='gray')
plt.title("Gradient")
plt.show()
strell = morpho.disk(6)
img_reco2 = morpho.reconstruction(morpho.opening(img_reco, strell), img_reco, method='dilation')
plt.imshow(img_reco2, cmap = "gray")
plt.title('Reconstruction from opening')
plt.show()
strell = morpho.disk(10)
img_reco = morpho.reconstruction(morpho.closing(img_grad, strell), img_grad, method='erosion')
plt.imshow(img_reco,cmap = "gray")
plt.title('Reconstruction from closing')
plt.show()
strell = morpho.disk(12)
img_erosion = morpho.erosion(img_reco2, strell)
plt.imshow(img_erosion, cmap='gray')
plt.title('Erosion')
plt.show()
strell = morpho.disk(20)
markers = peak_local_max(img_erosion, footprint=strell, min_distance=20, threshold_abs=3)
marker_img = np.zeros_like(img)
for k in range(len(markers)):
marker_img[markers[k][0]][[markers[k][1]]] = k
wat = seg.watershed(img_reco, markers=marker_img, watershed_line=True)
plt.imshow(wat, cmap=cmaprand)
plt.title("Watershed")
plt.show()